/**
 * @file CD.cpp
 * @author Lishijie (lsj1018845759@outlook.com)
 * @brief 利用Galerkin有限元法求解Preblem3.1.3，求解器利用Eigen中BiCGSTAB。
 *        计算结果在outflow layers（top）产生震荡，与Fig3.9一致。测试过P1和Q1。
 *  * @version 0.1
 * @date 2020-12-18
 * 
 * @copyright Copyright (c) 2020
 * 
 */
//g++ -o main CD.cpp -std=c++11 -I /usr/include/eigen3/ ./include/tinyxml2.cpp -g
#include "FEMSpace.h"
#include "Mesh.h"
#include "Element.h"
#include "Equation.h"
#include <Eigen/Sparse>
#include "Error.h"
#include "Eigen/IterativeLinearSolvers"
#include "MultigridSolver.h"
#include "Surface.h"
#include <cmath>
#define pi 4 * atan(1.0)
//========右端项==========//
Real f(Real* p)
{
    return 0;
}
//========外流===========//
Real wx(Real* p)
{
    Real x = p[0];
    Real y = p[1];
    return -sin(pi / 6.0);
}
//========外流==========//
Real wy(Real* p)
{
    Real x = p[0];
    Real y = p[1];
    return cos(pi / 6.0);
}
//========外流==========//
std::vector<Real> w(Real *p)
{
    return std::vector<Real>({wx(p),wy(p)});
}
//========bcf=========//
Real bc(Real *p)
{
    Real x = p[0];
    Real y = p[1];
    Real ans = 0.0;
    if(x == -1 || y == 1 ||(y == -1 && x <= 0))
        ans == 0.0;
    else if(x == 1 || (y == -1 && x > 0))
        ans = 1.0;
    else
    {
        std::exit(-1);
    }
    return ans;
}
//=======扩散系数========//
Real epslion()
{
    return 1.0/200.0;
}
typedef Eigen::SparseMatrix<Real> SpMat;
typedef Eigen::Triplet<Real> Tri;
typedef Eigen::VectorXd VectorXd;
int main(int argv,char *argc[])
{
    int n = 5;
    RectangleDomain * Domain= new RectangleDomain({{-1,-1},{1,-1},{1,1},{-1,1}});
    Mesh<2> * mesh = new Q1Mesh(Domain,{POW2(n),POW2(n)});
    Element<2> * element = new Quadrilateral_1_Element();
    int n_dofs_point = mesh ->n_dofs();
    int n_Dofs = element->n_Dofs();
    //=========组装stiff对流项扩散项==========//
    SpMat H(n_dofs_point,n_dofs_point);
    SpMat A(n_dofs_point,n_dofs_point);
    std::vector<Tri> TriList((mesh->n_element()) * n_Dofs * n_Dofs);
    std::vector<Tri> TriListM((mesh->n_element()) * n_Dofs * n_Dofs);
    std::vector<Tri>::iterator it = TriList.begin();
    std::vector<Tri>::iterator it_m = TriListM.begin();
    for (int k = 0; k < mesh -> n_element(); k++)
    {
        std::vector<int> Idx = mesh->NodeofEle(k);
        std::vector<Dofs<2> > temp(n_Dofs);
        for(int i = 0;i < n_Dofs; i++)
        {
	        temp[i] = mesh->DofsofIndex(Idx[i]);
        }
        element->SetDofsList(temp);
        for (int i = 1; i <= n_Dofs; i++)
	    for (int j = 1; j <= n_Dofs; j++)
	    {
            Real det = element -> det_Jacobi(0,0);
            Real a = 0;
            Real m = 0;
            for (int h = 0; h < element->n_GaussPnt(); h++)
            {
	            Real xi = element->GaussionPoint(h)[0];
	            Real eta = element->GaussionPoint(h)[1];
                Real p[2] = {xi,eta};
                a = a + element-> det_Jacobi(xi, eta) * element->GaussionWeight(h) * InnerProuduct(element->gradient(xi,eta,i),element->gradient(xi,eta,j));
                m = m + element-> det_Jacobi(xi, eta) * element ->GaussionWeight(h) * InnerProuduct(element ->gradient(xi,eta,j),w(p)) * element ->phi(xi,eta,i);
            }
		    *it = Tri(element->NdIdx(i), element->NdIdx(j), a);
            *it_m = Tri(element->NdIdx(i), element->NdIdx(j), m);
		    it++;
            it_m++;
	    }
    }
    A.setFromTriplets(TriList.begin(), TriList.end());
    A.makeCompressed();
    H.setFromTriplets(TriListM.begin(), TriListM.end());
    H.makeCompressed();
    A = epslion() * A + H;
    //===============组装rhs =======================//
    int n_GaussPnt = element->n_GaussPnt();
    VectorXd rhs = Eigen::MatrixXd::Zero(n_dofs_point,1);
    for (int k = 0; k < mesh -> n_element(); k++)
    {
        std::vector<int> Idx = mesh->NodeofEle(k);
        std::vector<Dofs<2>> temp(n_Dofs);
        for(int i = 0; i < n_Dofs; i++)
        {
	        temp[i] = mesh ->DofsofIndex(Idx[i]);
        }
        element->SetDofsList(temp);
        for(int i = 1; i <= n_Dofs; i++)
        {
            Real a = 0.0;
            for (int j = 0; j < n_GaussPnt; j++)
            {
		        double xi = element->GaussionPoint(j)[0];
		        double eta = element->GaussionPoint(j)[1];
		        double xi_ = element->Global_x(xi, eta);
		        double eta_ = element->Global_y(xi, eta);
                Dofs<2> temp({xi_,eta_});
		        a = a + element->det_Jacobi(xi, eta) * element->GaussionWeight(j) * element->phi(xi, eta, i)* f(*temp);
            }
            rhs[element->NdIdx(i)]+= a;
        }
    }
    //==============处理=========================//
    for(auto k:mesh->Boundary())
    {
        Dofs<2> bnd_point = mesh ->DofsofIndex(k);
        Real bnd_value = bc(*bnd_point);
        rhs[k] = bnd_value;
        A.coeffRef(k,k) = 1.0;
        for(Eigen::SparseMatrix<Real>::InnerIterator it(A,k);it;++it)
        {
            int row = it.row();
            if(row == k)
                continue;
            A.coeffRef(k,row) = 0.0;
        }
    } 
    //===============求解及可视化====================//
    Eigen::BiCGSTAB<SpMat> solver;
    solver.compute(A);
    VectorXd x = solver.solve(rhs);
    std::cout << " itertions is :" << solver.iterations() << std::endl;
    std::cout << " error is :" << solver.error() << std::endl;
    Surface<2> S(mesh,element,x);
    S.WriteVTKData("CD");
    S.WriteVTUData("CD.contour");
    delete mesh;
    delete element;
    delete Domain;
}